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PREFACE 


It is well known that a calculus-of —variations approach, to 
solving the Bolza form of trajectory optimization problems usually yields 
a nonlinear two— point boundary value problem in terms of the state and 

b 
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Lagrange multiplier variables. Closed— form solutions for problems of 
this type are difficult to obtain except for a few simple problems. As 
a result, recent work in trajectory '"optimization has focused on numerical 
procedures for obtaining solutions using high-speed digital computers. 
Particular interest has centered on a group known as the second-order 
methods . 

One such method is the Successive Sweep Kethod (SSM) . It uses 
the generalized Riccati transformation technique to bypass a direct nu- 
merical integration of the perturbation equations. The reasons such an 
approach has much potential appeal are presented in this study; however, 
because the SSM iterates on the control values over the interval of in- 
terest, considerable computer storage is necessary even for problems of 
small dimension. This storage is required to compute corrections to the 
assumed control. Furthermore, the Eulerian control is not obtained upon 
convergence . 

This research develops a new second-order numerical optimization 
-method, the Modified Sweep Method (MSM) . It requires very little computer 
storage and provides the Eulerian control. In addition, the properties 
and information contained in the Riccati transformation variables are 
preserved. 

Furthermore, this research also presents a new scheme for defi- 
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ning classes of numerical optimization methods. The Successive Sweep 
Method and the Modified Sweep Method are then discussed in terms of 

t 

differences arising because each falls into a different class. The Modi- 
fied Sweep Method is subsequently compared numerically to the Method of 
Perturbation Functions (OTP), both of which belong to the same class. 

The author extends thanks to Mr. Walt Williamson of The Univer- 
sity of Texas at Austin for many helpful discussions concerning the Apollo 
reentry problem, to Mr. I. J. Kim, of Lockheed Electronics Company, Hous- 
ton, for programming the plots, and both to Mr. Kim and Mr. Mike Frederick 
of The University of Texas at Austin for helping with the data. He also 
expresses gratitude to Dr. J. M. Lewallen of NASA/MSC and Professors A.M. 
Bedford and E.J. Prouse of The University of Texas at Austin for serving 
on the dissertation committee and helping with the manuscript. Special 
thanks is due to Mr. E. L. Davis of NASA/MSC for making this research 
possible, and for his personal interest and friendship which have pro- 
vided constant inspiration. 

The author expresses his indebtedness to Professor B. D. Tapley 
of The University of Texas at Austin, who suggested this research and 
provided many stimulating discussions while serving as committee super- 
visor. 

The author expresses deep appreciation to his family for their 
understanding and cooperation during the course of this research. 

Daniel Colunga 


iii 



A MODIFIED SWEEP METHOD FOR 


CONTROL OPTIMIZATION 
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Supervising Professor: B. D. Tapley 

A Modified Successive Sweep Method was devised which yields 
Eulerian solutions to two-point boundary value problems of control opti- 
mization. This was accomplished by requiring control satisfaction of 
local optimality over the entire time interval of interest while simul- 
taneously relaxing terminal transversality requirements on the Lagrange 
multipliers . 

The new method was tested successfully on several classes of 
problems including optimizing the roll program for an Apollo-type three- 
dimensional reentry ^trajectory so as to minimize a time integral of 
spacecraft heating and acceleration. 

This new method was-shown to require significantly less computer 
storage than the original Successive Sweep Method while requiring numeri- 
cal integration of fewer variables. In addition, the method was shown to 
possess rapid terminal convergence and a conjugate— point test capability. 
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CHAPTER 1 


INTRODUCTION 

Optimal control of spacecraft trajectories requires obtaining 
an optimal (maximum or minimum) value for an appropriate scalar quantity. 
This scalar quantity measures spacecraft performance and is called the 
perfoimance index. In addition, terminal conditions such as might be 
specified for intercept or rendezvous problems must often be satisfied 
simultaneously. 

After the problem has been formulated mathematically, several 

conceptual approaches are available to obtain the conditions required to 

$ 1 
. solve the optimization problem. Among the more usual approaches are the 

calculus-of-variations , dynamic programming, and Pontryagi^s Principle, 

The calculus-of- variations is considered here because specific 

optimal control problems can be considered as particular cases of the 

more generalized Bolza problem of the classical calculus-of-variations • 

j 

The powerful results associated with this classical theory thus are avail- 
able for attacking optimal control problems. 

The calculus-of-variations approach yields a nonlinear two-point 
boundary value problem for which closed-form solutions are usually not 
possible. Sophisticated numerical procedures have been developed because 
of the need to solve these problems. These procedures have become feas- 
ible in the last decade because of the development of large-scale digital 
computers . 

This chapter introduces the definitions and terminology used 
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throughout the dissertation. A brief history of the development of the 
numerical optimization methods is given also with interest centering on 
the second-order variational methods. The class of trajectory optimi- 
zation problems to be solved is stated, as well as the associated non- 
linear two-point boundary value problem obtained from using a calculus- 
of- variations approach. 

Chapter II discusses (1) the second-order variation approach 
used to solve the nonlinear two-point boundary value problem, (2) the 
general set of perturbation equations for second— order methods, and (3) 
techniques used to achieve an integration of these perturbation equations. 

Chapter III presents the Modified Sweep Method based on the 
generalized Riccati transformation. Chapter XV develops the linear feed- 
back control law for the Modified Sweep Method. The numerical results 
obtained using this new method are discussed in Chapter V, with conclusions 
and recommendations presented in Chapter VI. 

1.1 Definition of Terms Used 

The MSM (Modified Sweep Method) is obtained from the v SSM (Suc- 
cessive Sweep Method) by requiring that the control satisfy both local 
optimality and strengthened Legendre-Clebsch condition over the entire time 
interval of interest. This optimal control is then eliminated from the 
Hamiltonian for the problem and the restructured Hamiltonian used to obtain 

the nonlinear differential equations for both the state and Lagrange multi- 

1 

plier variables. As in the case of the SSM, the MSM then uses the 

generalized Riccati transformation to solve the linearized two-point 

* 

boundary value problem in terms of the state and Lagrange multiplier per- 
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tur'kations. This is done while simultaneously relaxing the terminal trans- 
vergality requirements on the time-dependent Lagrange multipliers. It is 
desirable to compare the proposed XSM to other existing numerical opti- 
mization methods. For this reason, a study was made of several well-known 
methods which appear in the literature. This author felt that these methods 
were representative of the properties contained in the set of variational 

oetfois for the numerical solution of optimization problems. It is empha- 

\ 

sizfd that only a representative portion of the total number of existing 
nunfrical methods has been used for this study. In addition , the second- 
ordf r methods intentionally have been selected more extensively than the 
first-order methods* The generalizations made, therefore, pertain only to 
tho^a methods contained in Section 1*2 on the historical development of 
numerical optimization methods* 

This study of the selected group of existing methods revealed a 
set of properties which can be employed to describe the characteristic 
f^^ures for each method* These properties have been used to specify, ar~ 
k^t.arlly, twelve classes of numerical optimization methods. Each method 
can then be identified as belonging to a particular class according to the 
f° -lowing properties: 

(1) the ORDER (first or second) of the theory upon which the 

method is based. 

(2) the APPROACH (direct or indirect) used By the method to 

compute the required corrections. 

(3) the ITERATION PROCESS (interval, boundary or hybrid) used 


by the method. 
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The following definitions are used: 

Definition 1: Order of the Method 

A method is described as first order if it is based only upon 
the theory of the first variation for a real functional- If a method is 
based upon the theory of the second variation for a real functional, it is 
described as a second order method. 

Definition 2: Approach of the Method 

A method is said to take a direct approach if the required 
corrections (state, control or Lagrange multipliers) are computed such 
that the performance index for the augmented variational problem is 
itself directly affected in some manner to expedite convergence. If a 
method chooses to compute the required corrections based on the set of 
first-order necessary conditions required for optimality with respect 
to the control, then the method is said to take an indirect approach. 

Definition 3: Iteration Process for the Method 

A method is said to use an interval value iteration process if 
the end result of a particular iteration is the computation of correc- 
tions to the variables (state, control, or Lagrange multipliers) over the 
entire time interval of interest. Methods which compute corrections to 
these same variables at a boundary only are said to use~a boundary value 
iteration process. If a method combines both an interval and boundary 
value iteration process, it is described as using a hybrid iteration pro- 
cess. 

Definition 4 : Convergence f or a Method 


(a) Control Function Iteration Method . Given an arbitrary 




numerical tolerance e, a control function iteration method is said to 
have achieved convergence if 

I Mh^II + l|H f || + n| | 1 6 \ 

where 

||H^|| = Max {Abs [H^(x,u,X,t) ] } t < t < t_ 

U u(t) U o t 

(b) Boundary Value Iteration Method , Given a numerical toler- 
ance e, a boundary value iteration method is said to have achieved con- 
vergence if 


1 IUJI 

+ 

I ImJ 1 

+ 



11 f 1 1 

f 


The symbols H, and are defined in Section 1.4* 

Using these definitions , Table I summarizes the twelve classes 
of numerical optimization methods extracted from the methods chosen as 
representative for the study. Reference numbers specify major studies in 
each class while the acronyms (SS M, etc.) identify the particular class 
for the three methods to be compared in detail. 



TABLE I 
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CLASSES OF NUMERICAL OPTIMIZATION METHODS 


CLASS 

ORDER 

APPROACH 

ITERATION PROCESS 

REFERENCES 

1 

1 

Direct 

Hybrid 

- 

2 

1 

Direct 

Interval value 

4,12 

3 

1 

Direct 

Boundary value 

- 

4 

1 

Indirect 

Hybrid 


5 

1 

Indirect 

Interval value 


6 

1 

Indirect 

Boundary value 


7 

2 

Direct 

Hybrid 


8 

2 

Direct 

Interval value 

19,26 

9 

2 

Direct 

Boundary value 

35 

10 

2 

Indirect 

1 

Hybrid 


11 

2 

Indirect 

Interval value 

23,27 (SSM) 

12 

2 

Indirect 

Boundary value 

10,11,16 

(MSM,MPF) 


As shown in Table I, the class of second-orddr methods which 
take a direct approach in computing the desired corrections and implement 
a hybrid iteration process was not represented in the methods chosen for 
the study. Furthermore, Classes 1,3, 4, 5 and 6, all first-order methods, 
were also not represented. This can be attributed to the fact that second- 
order methods were of primary interest in this study. 
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The Similarities and differences between the three numerical 
tlcn mcrhodr* (SSM, KSM and MPJ?) to be discussed are now obvious. 

1 are £^cond--order methods which use an indirect approach in 

mputi-^S the required corrections for the state, control or Lagrange 

variables- The SSM falls into Class 11 because it uses an 
cervai iteration process. Both the MSM and KPF fall into Class 12 
causes each uxes a boundary iteration process. 

* A hiiuorical development of the methods chosen for the study 
det£i>-vd nor for reference purposes - 

- H^g tcricE- Information 

Firsi-Order Methods - The first numerical procedure for solving 
~roi- c^rxmiiation problems which generated active interest was devel- 
d ifld^endarcly by Kelley* ^ and Bryson and JDenham** . Their research 
•ended *he except of steepest descent developed earlier by Courant^- 
Cla*>£ 2 meihod was based on the first-order variation of a scalar 
tio»£, wmih a control function assumed for the time interval of in- 

i 

.st- corrections to this control were then computed iteratively using 

rdiracy gradient technique. Applications showed that the method was 

to -implement and tended to convergence with even gross initial control 

\ 

sees*. 

The method, however, possesses two undesirable features. First, 
-■'ergeuce rate decreased asymptotically during the terminal stages 

Second, once convergence was achieved, the control obtained 
.. 'ri.tr.ia a numerical tolerance of the Eulerian control. 

Due to the first undesirable feature, numerical procedures to 
- the convergence rate flourished. These were all first-order 



methods and are not discussed in this research. Both features led to the 
development of the second-order methods which sought to increase conver- 
gence as well as provide the true Eulerian control. 

Second-Order Methods . Jurovics and McIntyre 11 solved the two- 

t 

point boundary value problem of trjaectory optimization by using the equa- 
tions which are adjoint to the linearized Euler-Lagrange equations. Their 
method was called the Adjoint Method (Method of Adjoint Functions), and 
extended the work of Goodman and Lance' 7 to allow for variable terminal 
time. The Adjoint Method is a Class 12 method wherein an indirect approach 
is used to compute the desired corrections while iterating on the initial 

boundary values of the Lagrange multipliers. 

Breakwell, Speyer, and Brvson 3 developed a "second variation" 

method (Class 10) to solve control optimization problems. Kelley, Kopp, 
and Meyer 13 also developed a "second variation" method similar to the 
previous one. Jazwinski 10 developed a modified adjoint method equivalent 
to the second-variation method of Breakwell, Speyer, and Bryson 3 by ex- 
tendiig the method of Jurovics and McIntyre 11 . Jazwinski *s method had 
the specific advantage, however, of requiring considerably less storage. 
Furthermore, it required less computer time in that fewer integrations of 

an equivalent set of equations were necessary. 

_McGill and Kenneth 20 developed the Generalized Newton-Raphson 

Operator Method for solving two-point boundary value problems . This method 
falls into Class 11, which uses the indirect approach linked with an 
internal value iteration process. A proof of quadratic convergence for 
the suchod had previously been given by these same authors 21 . This method's 
major drawback was the laborious manner in which corrections to the final 
time <alue were computed. 



An alternate approach, the Modified Generalized Newton-Raphscn 
cchod-j using the. same method V7a$ developed by Long 17 * To eliminate the 
.xwa fd handling of free final time, a change of independent variable was 
sed £k>r the frea-f inal-txme corrections. Another method based on the 
swtort-Raphson approach, but incorporating a' better technique for compu- 
ting the free-final-time corrections , is the Modified Quasilinearisation 
echod developed by Lewallen 16 . Sylvester and Meyer 36 have also used the 
awtoh^Rephson approach, calling it quasilinearization. 

A method based on the theory of both the first and second varx- 
.txond was devised by Merriam 26 . This is a Class 8 method in which a 
irect approacr. is taken for computing corrections to the control functions 
ssumed throughout the time interval of interest. This particular method 
as instrumental in the development of the successive sweep method discus- 
ad in the next paragraph. 

McReynolds and Bryson 26 introduced the successive sweep method 
•r solving optimal control problems. Although the method is a second- 
rd»:r method, the Eulerian control requirement is relaxed. The method is 
-ed on the generalized Riccati transformation and falls into Class 11 
ible X) . A similar method called successive approximation was developed 
.‘fitter 27 , ha, also showed the formal equivalence of this method to 
-on’s Metrics . , 

Lewallen^ 5 also introduced the Method of Perturbation Functions 
*r), based or previous work by Breakwell, et al . , 3 and Jazwinski 16 . The 
.nod falls into Class 12. Lastman^ 1 * has shown the equivalence of all 
•se methods to Newton's Method. ' / 

Sutherland and Bohn 35 have recently developed a method which falls 
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into Cl^ss 9 , which, uses a direct approach to compute boundary corrections 
for the initial values of the multipliers. 

Mayne 19 developed a second-order method which falls into Class 8 
(Table I). However, a dynamic programming technique is used to attack the 
optimization problem instead of a calculus-of-variations technique. 

Jacobson 9 extended the Class 11 features into a new second-order 
algorithm through use of a differential dynamic programming technique. His 
method generalizes the successive sweep methods of McReynolds 29 and 
Hitter 27 . 

I 

The development of the MSM completes the historical development 
for the numerical optimization methods chosen. 

As was mentioned previously, it is now desirable to compare the 
MSM to both the SSM and the MPF. Toward this end, a general class of con- 
trol optimization problems is first chosen. This class of problems is 
presented in the next section. 


1.3 Class of Control Optimization Problems to be Solved 

Posed as a special form of the Bolza problem from the calculus- 
of-variations (see Bliss 1 ), the general class of control optimization prob- 
lems to be solved is stated as follows: 

In the time interval t < t < t. , find an m-vector of control 
variables u(t) to minimize the real functional. 



t 

o 


J(u) 


G(x f ,tf) + 


Q(x,u,t) dt 


( 1 ) 
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subject to the n-vector of differential constraints, 


x - f(x,u,t) 


while satisfying the p— vector of known initial conditions 


N(x ,t ) - 0 

o o 


and the q— vector of desired terminal conditions 


K(x f ,t f ) = 0 


( 2 ) 


(3) 


(4) 


The control and state variables in the following discussion are 
assumed to be defined on completely open regions and thus are not subject 
to inequality constraints. 

1.4 Associated Nonlinear Two-Point Boundary Value Problem 

Proceeding in the usual calculus-of-variations manner for solving^ 
the Bolza problem of control optimization stated in the previous section, an 
augmented functional denoted as I is first formed. This augmented func- 
tional has the property of being formally equivalent to the original func- 
tional; it incorporates the desired auxiliary conditions through the use 
of Lagrange multipliers. To form this augmented functional, the n-vector of 
Lagrange multipliers A(t) and the p and q vectors of constant Lagrange 
multipliers y and v adjoin the desxred auxiliary conditions to the ori- 
ginal functional as follows: 


I = J + p T N + V T M + 


X (f - x) dt 


(5) 
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For convenience of notation , this functional is rewritten as 


where 


I = P + V + 




(H - X x) dt 


( 6 ) 


P = P(x f ,v,t f ) = G(x fJ t f ) + v^Cx^t.) 


f’ f' 


V = ^(x o> y,t o ) & A(x o ,t o ) 


H = H(x,u,X,t) = Q(x,u,t) + X f(x,u,t) 


and 


x = x(t) 


u = u(t) 


X = X(t) 


The scalar H is the variational Hamiltonian for this class of problems. 


Necessary Conditions . The set of first-order necessary condi- 
tions which must be satisfied by the extremal control for the augmented 
functional of the type above is obtained by requiring the first variation 
of this functional to vanish. These conditions are well documented in the 
literature (for example Bliss 1 , Hestenes 8 , Pontryagin 30 , et al . , Tapley and 
Lewallen 38 ). In summary, these conditions are 


x - H^(x,u,X,t) = 0 


t Q < t < ■{ X + H^(x,u,X,t) - 0 


H (x,u,X,t) = 

u 


( 7 ) 

( 8 ) 


0 


(9) 
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* - C o < 


i? = (T + X T ) = 0 

o x 't 

o 


N(x ,t ) = 0 

o’ o 


(105 

0 - 1 } 


"o S <V H >t - 0 

o 


( 12 } 


t = t. 


ll * (P - x T ) 

f x y t f = 0 


M(x f ,t f ) = 0 


(13) 

0 - 4 ) 


fi f ■ < p t + H h t “ 0 

- i— 


(15) 


v 

Equations (7) through (9) constitute 2n+m Euler-Lagrange equations for 
this class of problems. Equations (11) and (14) are the p+q specified 
initial and final values of the problem state variables. The remaining 
equations form the 2n+2 set of classical transversality conditions from 
the calculus— of- variations. 

The control optimization problem thus is posed as a nonlinear 
two-point boundary value problem for the 2n-hn variables x(t), u(t) , 
and X(t) and the p+q+2 parameters y, v, t Q , and t^ in terms of 
the 2n differential equations (7) and (8) , the m algebraic equations 
(9), and the 2n+p+q+2 conditions in equations (10) through (15). 
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It is assumed that the initial time t and n values ox the 

0 

initial state x(t Q ) = x q are specified. Equations (10) and (12) are 
then identically satisfied and therefore disregarded in subsequent dis- 
cussions. 

1 

Sufficiency Conditions . To ensure that the control satisfying 
these first-order necessary conditions does indeed generate a weak mini- 
mizing solution, the second-order variation of the augmented functional 
I(u) must be positive everywhere in the interval of interest when it is 
evaluated along an extremal trajectory (Gelfand and Fomin 6 ). This re- 
quirement leads to the following additional set of second-order conditions 

1. Strengthened Legendre-Clebsch Condition 

The strengthened Legendre-Clebsch condition must be satis- 
fied everywhere in the interval of interest. Specifically, 
for any arbitrary change in control <5u(t) , 

Su T H 6u > 0 (16) 

uu 

is required. 

2. Jacobi (Mayer) Conjugate Point Condition 

The Jacobi (Mayer) conjugate point condition must be satis- 
fied everywhere in the interval of interest'. This requires 
that no two points exist in the interval t < t < t f which 
are conjugate to one another. 

The following restrictions on the definitions presented in Sec- 
tion I are subsequently assumed m this research: interval iteration 

corresponds to control function iteration and boundary iteration corres- 
ponds to Lagrange multiplier iteration. 


CHAPTER 2 


THE PERTURBATION EQUATIONS 

The second-order variational methods seek to solve the nonlinear 

$ 

two-point boundary value problem associated with trajectory optimization 
by solving an equivalent linearized problem in terms of perturbations in 
the problem variables. The six classes of second-order methods are dis- 
tinguished by approach and iteration process. Furthermore, each method 
for a given class is distinguished by the technique used to perform the 
numerical integration of the perturbation equations which are obtained 
from a first-order perturbation of the Euier-Lagrange equations, viz.. 
Equations (7) through (9). These perturbation equations can assume one of 
two forms with each form patterned by the iteration process selected for a 
given method. 

As in the case of the SSM (Successive Sweep Method) , if a con- 
trol function iteration process is used, a "PE" scheme is used to obtain 
the perturbation equations. The acronym "PE" designates the following 
procedure: Perturb the Euier-Lagrange equations and then Eliminate the 

control perturbations . If a Lagrange multiplier boundary value iteration 
process is used, an "EP" scheme is used to obtain the perturbation equa- 
tions. The acronym "EP" similarly designates the following procedure: 
Eliminate the optimal control and then Perturb the resulting Euier-Lagrange 
equations. These two schemes for generating the perturbation equations 
for second-order methods are detailed below. 

2.1 The PE Scheme 

The PE Scheme is associated with control function iteration 


15 
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schemes such as the SSM, It proceeds by first perturbing the first- 

order necessary conditions for stationary control. Then the perturbations 

in the control are elininated from the perturbed Euler-Lagrange equations 

for the state and Lagrange multipliers and the appropriate transversality 

conditions. This scheme tacitly assumes that the matrix tt is nonsin— 

uu 

gular everywhere in the time interval of interest. The set of perturba- 
tion equations associated with this scheme is outlined in Appendix A and 
are summarized here as 


where 
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A = 


-H, H _1 H 

AU UU UX 


+ H 


Xx 


B = 


-H, H _1 H , 

AU UU UA 


C = -H H~ 1 H + H 

XU UU UX XX 


V = 
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H. H 6H 
Xu uu u 


w - 


-I T 
H H <$H 

XU UU U 


and H = K(x,X,u,t) while the Hamiltonian partials H , etc., are 

evaluated on the known trajectory. The known trajectory is called a nomi- 
nal (or reference) trajectory. 

Note that for this scheme, the system of perturbation equations 
assumes the form of a set of inhomogeneous first-order linear differential 
equations with time-dependent coefficient matrices. 
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2.2 The EP Scheme 

The EP Scheme is used by the boundary value iteration methods. 

T 

It proceeds by first using the control optimality condition H = 0 and 

T 

the strengthened Legendre-Clebsch condition 6u H 6u > 0 to eliminate 
the control from the Euler-Lagrange equations for the, state and Lagrange 
multipliers, as V7ell as the appropriate transversality conditions. The 
revised set of equations are perturbed then to obtain the following homo- 
geneous set of first-order linear differential equations. 
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(IS) 


where 



and u(x,X,t) is the control obtained using the conditions that 0 

and H >0. It is imoortant to note that this scheme involves the as- 
uu 

sumption that the Hamiltonian H for the particular problem is structured 

such that H T - 0 and H >0 can be used to obtain the explicit rela- 
u uu 

T 

tion u(x,X,t). Implicitly, such a relation is assured if *= 0 and 

H > 0 ; however 5 such an explicit relation may be impossible to obtain 
uu 


for sorse problems* 



2.3 Integrating the Set of Perturbation Equations 


It has been pointed out previously that second-order methods 

t 

within a given class differ only in the technique that is used to inte- 
grate the set of perturbation equations. The two techniques presently 
available for accomplishing this integrating are detailed below. 

Explicit Integration . Methods using this technique choose to 
integrate directly the perturbation equations to obtain the perturbed 
values over the interval of interest. Some investigators have found (see 
for example the work of Merriam 26 ) that these methods suffer from numer- 
ical instabilities. Instability here is used in the sense that small 
errors in numerical precision will become exponentially very large over a 
long interval of numerical integration. The nature of these instabilities 
is associated with the numerical integration for coupled systems of linear 
differential equations which have split boundary conditions and admit both 
increasing and decreasing exponential solutions. 

Implicit Integration . Methods which presently use this technique 
are based upon a transformation process such as the generalized Riccati 
transformation. The technique consists of bypassing direct integration of 
the perturbation equations, and integrating a set of auxiliary variables. 
These in turn can be used to compute the perturbed values for the variables 
in the perturbation equations. Several advantages are claimed for this 
technique. First, the differential equations for the new auxiliary vari- 
ables have been reported to be more stable numerically than the original 
perturbation equations. Second, these auxiliary variables contain addi- 
tional intrinsic information about the optimal trajectory for the problem 
being solved. 



A strong case concerning increased numerical stability in first- 

♦ 

order linear two-point boundary value problems has been made by Rybicki 
and Usher 31 . However, work by Williamson 39 and this author has revealed 
that problems which have large differences in sign and magnitude for the 
eigenvalues of the coefficient matrix in the linear system of equations do 
not behave well numerically with the generalized Riccati transformation 
technique. This author’s opinion is that a valid generalized statement is 
yet to be made concerning the numerical stability properties of the Riccatf. 
transformation technique. 

That the auxiliary variables could contain additional intrinsic 
information certainly proves to be true. The earlier methods lacked in 
.one respect: after convergence had been achieved, they required that post- 

convergence procedures be used to test the Legendre-Clebsch condition and/ 
or the Jacob i-Mayer conjugate point condition. Those methods which used 
the "EP" Scheme to generate the set of perturbation equations automatically 
took the Legendre-Clebsch condition into account when eliminating the con- 
trol from the set of Euler-Lagrange equations. However, the Jacobi-Mayer 
conjugate point condition still must be tested. This condition was often 
ignored and the converged solution was assumed to be a local optimum. 

However, using the generalized Riccati transformation technique 
on the perturbation equations provides the advantage of additional infor- 
mation for the current reference trajectory. Information is contained 
among these auxiliary Riccati variables for testing the Jacob i-Mayer con- 
jugate point and abnormality conditions from the calculus— of— variations 
(McReynolds 23 ) . It is well-known that the existence of a conjugate point 
precludes a trajectory from being optimal. The existence of such a con- ' 
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jugate point can be automatically detected continuously during the back-* 
ward sweep process. This is accomplished by use of the fact that the 
matrix solution to the Riccati differential equation becomes unbounded at 
a conjugate point. ^ 

I 

On the other hand, the abnormality condition is equivalent to a 
certain matrix of these auxiliary Riccati transformation variables be- 
coming singular at the initial time. This information is important be- 
cause such a condition is tantamount to the inability in making correct- 
ions for values of the terminal constraints. This abnormality condition 
occurs for the Bolza problem if the boundary conditions at the final time 
are not linearly independent (McReynolds 23 ) . 

This leads to speculation concerning additional information 
about the reference trajectory which might be contained in the other Ric- 
cati transformation variables, either individually or in some combined 
form. To this author’s knowledge, little work has been done in attempting 
to extract such additional information. 

2.4 The Generalized Riccati Transformation Technique 

The generalized Riccati transformation is a transformation which 
changes the original two-point boundary value problem in terms of the 
coupled linear system of differential equations to an initial-value prob- 
lem having uncoupled variables and boundary conditions. This initial 
value problem is now stated in terms of -n original problem variables, 
and in the general case, a total of |n(n+i) + q(q+l)j /2 + j\nxq) + 

3(n+q) 4 - 2j auxiliary Riccati variables where n is the number of state 
variables and -q is the number of terminal constraint. Since the coupled 
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system of perturbation equations is integrated implicitly by integrating 
these auxiliary variables, it was expected that the differential equations 
for the auxiliary variables would be numerically more stable than the ori- 
ginal equations. Furthermore, Breakwell and Ho^ have shown agreement with 
McReynolds 23 m that the conjugate point condition is related directly to 
the boundedness of an (nxn) matrix of Riccati variables which must satisfy 
a matrix version of the scalar Riccati equation over the time interval of 
interest. 

This transformation approach proceeds to solve the original non- 
linear two-point boundary value control optimization problem in the fol- 
lowing manner. A solution to the original nonlinear problem is assumed, 
and the corresponding terminal conditions are obtained. In general, these 
conditions are not satisfied to within the specified error tolerances. 
Desired changes, in these terminal conditions are specified, and the gen- 
eralized Kiccati transformation is used then to generate a linearized field 
of solutions about this assumed solution. The transformation allows the 
specified changes in the set of terminal conditions to be mapped back to 
the initial time, when the particular member of the field that also satis- 
fies the initial conditions is selected. A new solution to the original 
nonlinear two-point boundary value problem is then computed using the lin- 
earized corrections, which are obtained through use of the auxiliary Ric- 
cati variables. As before, the new solution does not satisfy the desired 
terminal conditions exactly due to the linearity assumptions. However, the 
process can be applied iteratively until the desired terminal conditions 
are satisfied to within a suitable error tolerance. 

Historical Background. A Riccati transformation technique was 
first used by Gelfand and Fomin 6 in their successive sweep procedure of 
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solving two-point boundary value problems for linear inhomo geneoug^ 
of second-order differential equations, The same transformation^: 
eralized and discussed for systems of first-order equations by Ryhiii 
Usher* 32 McReynolds 23 * 2tf and Hitter 27 used the successive 
with the generalized Riccati transformation technique to solve' trie * 
linear two-point boundary value problem of control optimization* -gSc* 
and Lee 32 have developed a Newton-Raphson method which uses the Ric^ ~ 
transformation technique. Speyer and Byrson 34 have extended the c^gg* 
of the Riccati variables for the case when some of these variables^&''C^ 
unbounded. Narha and Berry 28 and Omicioli 2 ^ have applied the 
sweep method of McReynolds to the shaping of optimal finite-thsg^^g^ 
transfer trajectories for which the control function is characftep^gy 
discontinuities. McGregor 22 has used the same method but has 
modifications to handle problems with inequality constraints 
the control explicitly. Most recently, Longmuir and Bohn 1 8 have 
how this technique can be used with any second-order method. 

Analytical Development . The generalized Riccati transxiz^Ssr 
for the linearized control optimization problem can be written i 
‘form as 




SxCt)’ 

dM f 

- R(t) 

dv 

d^ £ 


dt- 
» , 


+ p(t) 


( 19 ) 



2B 

where dM^, dft^, dv> and dt^ are constants for a particular iteration 
and 



T 

where K, E, y respectively map given state perturbations 6x(t) into 

changes in the multipliers SA(t) > terminal state dissatisfactions * 

I 5 

and terminal Hamiltonian transversality dissatisfaction and 

f 

K(t) is an nxn symmetric matrix 

E(t) is a qxn matrix 

y(t) is an n-vector 

T 

Also a F, z y respectively map changes in the multipliers dv into 

changes in the multipliers 6A(t) 5 terminal state dissatisfactions dM„% 

* r 

and Hamiltonian terminal dissatisfaction d$J- s and 

f 

D(t) is an nxq matrix 

F(t) is a qxq symmetric matrix 

z(t) is a q-vector 

The scalars SL, g and s respectively map changes m the final 
time dt f into changes in the multipliers 6\( t) , terminal state dissat- 
isfactions dM f and the dissatisfaction in the terminal Hamiltonian, da . 

’ l 
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The quantities n, £ and 4 respectively map changes in termi— 

4 

nal local optimality dissatisfaction or terminal transversality dissatis- 
factions by the multipliers into changes <5X(t), dM„, and dSK . 

Differentiating the generalized Riccati transformation (Eq. 19) 
with respect to time gives 

6x <5x 

R dv + R 0 + p (20) 

dt f 0 

k \ d 

Expanding the middle term on the right and transposing to the left yields 

I -K 

n 

0 -E 

0 -y 

Using the perturbation equations. Equations (17), to eliminate 61 and 

• 1 

6x leads 'to the following expression 





( 22 ) 
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Using the first row of the Riccati transformation, Equation (19), the 
following relation can be readily obtained. 
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Substituting Equations (23) and (24) back into Equation (22) gives 
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Multiplying and collecting terras for 
the following equations must hold 


where 


and 



Performing the matrix multiplications 
tions for the Riccati variables. 
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From Equation (28) , the following rates of change for the Ric- 
cati variables are found to be equal: E^(t) = D(t), y(t) = 2,(t) and 

z(t) = g(t). If, then, at the terminal boundary E (t^) = D(t^), y(t^) - 
f(t^) and z(t^) = g(t^), the following will be true: E (t) = D(t), 
y(t) — 2.(t) and 2 (t) = g(t). This means not only that the matrix of 
Riccati variables R(t) given in Equation (19) is symmetric but also that 
Equation (26) itself is also symmetric. In this case, 


R = -SWS 


where S and W are defined on page 26. 


( 30 ) 


Terminal Boundary Conditions . The derivations of the general- 
ized set of terminal boundary conditions for the Riccati variables are 
presented in Appendix B. In summary, these boundary conditions are ob- 
tained from 
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XntepygZ Value Process , Methods using this process start V 

assuming a control function over the time interval of interest* In 
T 

general then, H^(t) ^ 0* For such methods, the Lagrange multiplier * 
can he made to satisfy the transversality conditions identically; 

Z^ ~ 0. The Successive Sweep Method is an example of a contrqi,Jiu[£fa. 
iteration process* The terminal boundary conditions for the Jg 
variables are then obtained from 


and 


where 


“f 

r 

dM, 

_ 

f 


dfl r 



V 


CP ) 
xx f 

M 


T . 

0t £ + T 
f 2 



• T 

M f 


+ T 

3 



t = f-H H _1 (H + H .P )> 

2 C u uu ux uX xx 

t = [— H H -1 H .M T 1 

3 ( u uu uX xj 

f 

T » f-H H'" 1 H .a,’) 

( u uu uX fj t f 

t - fn k~ ;l <sh t 'I 

5 [uuu ujt f 


and 



m 


NoteVSrat the terminal matrix R(t,.) of Riccati variables is 

r < 

symmetric only when the control function satisfies the local optimality 
T 

condition that H (t- ) = 0. This same conclusion has been drawn by 
u f J 

McGregor 22 . However, this contradicts the results presented by McReynolds 
and Bryson 25 and Mitter 27 . This argument needs to be resolved. 

T 

Boundary Value Process . These methods assumed that H^(t) = 0 , 
where £^ft *'-»-Jn this case, initial values are guessed for the 


Lagrange mui'ti^sers, and subsequent corrections are computed iteratively. 


In general the^,_ E^ ^ 0: likewise, dE_^ 4 0- The terminal boundary 

conditions in this case are obtained as 



(35) 


Note then that there is a basic difference in philosophy between a control 

function and a boundary value iteration process. In the first case, the 

1 

optimality condition that = 0 is relaxed at each point in the inter- 

val to ensure satisfaction of the transversality condition E^ = 0 by 

the Lagrange multipliers. In the second case, the optimality condition 
T 

= 0 is satisfied at each point in the interval while the terminal 
transversality E^ = 0 is relaxed on the multipliers. 
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THE MODIFIED SWEEP METHOD 

Merriam 26 and Mitter 27 have pointed out that bound ary- condition 
iteration methods have certain programming advantages ; viz . , computer 
logic is relatively simple, and programming storage requirements are 
small. Furthermore, accurate trajectories are obtained in problems where 
these methods are successful. Experience has shown that such methods have 
one main disadvantage, vi 2 . , the numerical instability mentioned previously. 
The nature of this instability has been discussed by several researchers, 
among them Rybicki and Usher. 31 Since the Riccati transformation tech- 
nique attempts to circumvent this problem by dealing with new uncoupled 
variables, this approach enhances the desirable features already known 
about boundary iteration schemes. , 

3.1 Differential Equations 

~ u 

For the modified sweep method, it is assumed that the Euler- 

Lagrange equations are satisfied over the time interval of interest* Fur- 

T 

thermore, the optimality condition that = 0 is assumed to yield an 
explicit expression for the control m terms of the other variables. The 
Legendre- Cl ebsch condition is then used to yield the extremal control 

& 

u = u(x,A,t) ’ - (36) 

This expression can now be used to eliminate the control from the original 
nonlinear Euler-Lagrange equations for x and X as well as from the 
appropriate trans vers air ty conditions. The set of first-order necessary 
conditions can now be rewritten as 
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t < t < t c 
o - - f 


t = t 
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t = t. 


• ~T 

x - H^(x,A,t) = 0 

(37) 

X + H^(x, A ,t) = 0 

(38) 

N(x o ,t o > = 0 

(39) 

T ' T 

z f - (P x- X) ■ 0 

(40) 

M(x f ,t f ) = 0 

(41) 

= (P fc + H) = 0 

(42) 


Equations (37) and ( 38 ) are 2n equations for the 2n unknowns x(t) 
and A (t) . Equations ( 39 ) through (42) are 2n+q+l conditions for the 
2n unknowns x(t) , A(t), and the q +1 unknown parameters v and t^. 
These equations constitute the familiar nonlinear two-point boundary value 
problem. A first-order perturbation of the nonlinear Euler-Lagrange equa- 
tions is now considered. This yields the following homogeneous linear 
system of equations (see Section 2 . 2 ). 
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As was mentioned on page 27 , the differential equation for the 

n+q+1 matrix of Riccati variables, R(t) , will be symmetric if the terminal 

T 

boundary values are such that R(tp = R (t^). In the next section it 
is shown that R(t^) will always be symmetric for the MSM. This 
reason, along with the fact.. that an EP scheme is used by the MSM to obtain 
the perturbation equations given in Equation (18) , give the following dif- 
ferential equations for the MSM Riccati variables: 
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Boundary Conditions 


Boundary conditions for these equations are obtained by a first- 
order perturbation of equations (39) and (40) 


6x(t o ) = 6 x q (46) 

+ M x f dv + “f dt f - • (47) 

Equations (43) represent 2n equations for the 2n unknowns 6x(t) and 
5X(t). Furthermore, equations (46) and (47) give 2n split-boundary con- 
ditions for these variables in terms of the 2n known parameters 
and 61^ plus the q+1 additional unknown parameters dv and dt^. The 
required additional q+1 conditions are obtained by also performing a 
first— order perturbation of equations (41) and (42) . It is shown in 
Appendix B that this procedure yields 
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In matrix form, these boundary conditions can be summarized as 


ex; 


'(P ) 

XX f 

T 

M 

X f 



r *> 

6x f 


r \ 

- d£ f 

dM. 

__ 

M 

0 

M. 


dv 

4 - 

0 

f 


X f 


f 





dfl f 


T 

•T 

M 

n f 

(B f + x, ) 
t 4 J 


dt f 


T 

5 J 


( 50 ) 


Note that this coefficient matrix is symmetric* Furthermore, if n values 
of the state are specified at t^, = dE^. - 0 . For this case, the ter- 

minal boundary conditions reduce to the homogeneous form 


ffixJ 

t 

X 


dM. 


f 


dft. 


f 

J 

L 


XX £ x f 


M 


a. 


0 M, 


• T 

M f 8 f 


> 

(x ^ 

Sx_ 


r 


dv 


dt_ 

J 



( 51 ) 



as 

To summarise. Equations (43), (46), and (50) constitute the linear fir, 
order, two-point boundary value problem for the 2n functions <Sx(t) , 
6A(t), and the q-KL parameters dv and dt, in terms of the 2n+q+l 
specifiable parameters , SA^ , dM^ , and dft^. The c 

procedure followed by the MSM will now be outlined. 

3 * 3 Computational Algorithm 

The modified sweep method can be implemented as follows: 

Assume n-fq-KL values for XCt^) s V and t^ . 

Integrate the nonlinear Euler-Lagrange Equations (37) anS^f£J 

forward from t to viz- , 

o V , 

x » H^(x, A, t) 

* = -H3(x,X,t) 

A 

Step 3 - Test the error norm 

| [Error! | = i ! £ f II + I t M f I 1 + « f f e - 

If this criterion is satisfied, exit; otherwise, continue to 


Step 1 — 

Steo 2 — 


the next step. 
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Step 4 


Step 5 


Step 6 


Set the terminal boundary conditions for the Riccati vari- 
ables 
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- Integrate the Riccati variables backward from the final to 
the initial time using the differential equations 
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Compute the n-fq+l corrections 6A q , dv anci dt^ using 
the generalized Riccati transformation (Equation 19), eval- 
uated at the initial time. These corrections are 
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Step 7 


dv - (F - gs" i g T )^ 1 ( (dM f - 5 ) -gs 1 (d« f - *) 

o 

- (D T - gs _1 *, T ) Sx )t 

o 

dt » s“ 1 ( (dfi, - <J>) -g T dv -£ T 6 x ) t 

t t o t ° 

and 


6X « f K6x + Ddv + A dt- + n L 

O K I J L 

where it has been assumed that 
dlC « -t M* 

f t 

* 

dft^ = -e £2^ 

* 

0 < e < 1 

Repeat from Step 2 using the new values 


x i+ 1 (t o ) - 

xX) 

+ 

i+l 

V = 

v 1 + 

dv 

r ■ 

fc f + 

dt 


o 


3.4 Computational Advantages 

The computational advantages of the MSM over the SSM are 
that it has to integrate n+q less variables and requires considerably 
less storage than the SSM* The exact comparisons are shown in Tables 
II and III. Table II shows the number of variables which must be inte- 
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grated by each method in forward and backward directions. Table III 
shows the storage requirements for each method. Only 3n + 2(q+l) quan- 
titles need to be stored in the MSM case. For the SSM, however, 

M (n (q+3) -rnvfn (n-KL) / 2) quantities have to be stored where M is the total 
number of points in the integration interval. A quick check ior a 
typical reentry problem with M = 1,000, n = 6, q-3, and m - 1 
shows the MSM requires storage of 26 quantities while the SSM must 
store 58,000 quantities. 

It is speculated that use of the SSM for large complex prob- 
lems such as the Apollo 3-D reentry will require fixed step— size integra- 
tion routines with a large enough step” size to remain within the computer 
storage limitations. Furthermore, the large step size may lead to unsat- 
isfactory numerical accuracy. 
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TABLE II ' ? 

MSM AND SSM VARIABLES TO BE NUMERICALLY INTEGRATED 


STANDARD SWEEP METHOD 

MODIFIED SWEEP METHOD 

Forward : 

x - n 

Forward: x - n 

Backward: 

X - n 

1 

X - n 


n(n-KL) 
K 2 

Backward: K - ^ 

i ^ 


D - nxq 

D - ixxq 


& ~ n 

l - n 


ri - n 

n - n 


F _ 4<9±Ii 

F _ qCq±i). 


2 

2 


, g - q 

g - q 


s - q 

c - q 


s - 1 

s - 1 


rH 

1 

-e- 

<p - 1 


y - n 



z - q 


Totals : 

• 

Totals : 

n (q+5 ) + 

C3q+2) 

n(q+4) + (2q+2) 


(n(n-KL) + q(q+l)) 

+ (n(n+l) r q(q+l)) 



Difference: (n+q) less variables 



Note: If all values of the term- 
inal state are constrained, i.e., 
if q = n, Lhen n = 0, £ = 0, 

4> = 0 and the difference increases 
to 2(n+q)+l less variables. 




TABLE III 

COMPUTER STORAGE REQUIREMENTS 


40 


STANDARD SWEEP METHOD 

MODIFIED SWEEP METHOD 

At every point in the integration 

At the initial and final points only 

interval, the following values 

the following values must be stored: 

must be stored: 


x(t) : n 

x - n 
o 

u(t) : m 

A - n 
o 

X(t) : n 

V - q 

K(t) . n< " +1) 

t f - 1 

D(t) : nxq 

Z £ - n 
f 

£(t) : n 

M f - q 


5 f -‘ 1 " 

Let M = total number of points in 

Only 3n+2(q+l) quantities need to 

the integration, then 

be stored from iteration to iter- 

M ^n(q+3) + m + 

ation. 

■ Compare 26 quantities with 58,000 

quantities have to be stored. A 
typical reentry problem has 

for the Apollo reentry problem. 

M = 0rder(l,000) , n = 6, q = 3, 


and m = 1. 


Thus, 58,000 quantities must be 
stored over the integration 
interval. 



CHAPTER 4 


THE MODIFIED SWEEP METHOD GUIDANCE SCHEME 

Initial conditions for dynamical processes are difficult to 
control in actual problems. Errors often occur which may be due to in- 
ternal mechanical causes such as premature cutoff by a thrusting rocket 
motor. Regardless of where these errors occur, they have the cumulative 
effect of causing a deviation from the intended optimal path. It is then 
desirable to use the known information about the path to recompute a new 
control program to accomplish the mission objectives . This is done by 
determining the control function corrections <5u as a function of the 
state perturbation; i.e., <$u = 6u(6x,t). This is a guidance problem in 

optimization theory. The guidance relations are now derived using the 
MSM. 

The MSM assumes that from the local optimality condition 

T T 

= 0 and the strengthened Legendre-Clebsch condition 6u H 6u > 0 , 

the minimizing control u(t) can be expressed as an explicit function of 

the state and Lagrange multiplier variables; i.e., 

u = U(x,l,t) (54) 

Perturbing this expression to first-order gives 

6u = U x 6x * (55) 

The generalized Riccati transformation is 


41 
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<5A = K<$x + 

D dv 

+ 

A dt 

f 

+ 

D 

(56) 

dM f = D T 6x + 

F dv 

4~ 

g dt f 

+ 

5 

(57) 

dft^ = A^6x + 

T 

g dv 

4- 

s dt £ 

4- 


(58) 


Solving the last too equations simultaneously yields 


where 


and 


dv = 

IT $X + 

ir y 4 ir z 


11 

12 

13 

dt. * 

iT 4x i 

it y -f ir z 

f 

21 

22 

23 


y = (dM f - S) 

z = (dS2 f - <*>) 

A = (F - gs g ) 

A A/ -l.T _T S 
ir = A(gs A - D ) 

11 


12 



r 

13 

A 

, -1 
-Ags 

r 

21 

A 

-!/ T 

-s (g IT + 

11 

r 

22 

A 

to 

-1 T 
-s g it 

12 

T 

23 

A 

-s" 1 (g 1 ^ 3 -l) 


(59) 

(60) 


"Using Equations (59) and (60) to eliminate dv and dt f from Equation 
(56) gives 

6 A = (K + Du + Air )6x + (Dir + Air )y 

11 21 12 22 


+ (Dir + Ait ) z 
13 23 


( 61 ) 
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Substituting this expression into Equation (55) then gives 

<5u - fll + U. (K + Dir + Jin )}<$x 

X A 11 21 

+ U. ((Dn + Jlir )y + (Bn + Jin )z} (62) 

12 22 13 23 v ' 

This equation represents the linear-feedback control law for continuously 
correcting the optimal control program for a given perturbation Sx(t) in 
the vehicle state. The new control program is then obtained by adding 
these corrections to the converged control program. Note that for a 
given perturbation <5x(t) , the coefficients are obtained by integrating 
the Riccati variables forward using the initial values corresponding to 
the converged solution. 

In guidance work it is desirable to know the control corrections 
in terms of a known time-dependent matrix and the initial vehicle state 
perturbations, i.e. , , 

<5u(t) = L(.t) <5x(_t o ) ( 63 ) 

where L(t) is an explicit relationship between the Riccati variables 
from the converged optimal trajectory. Attempts to yield a relation such 
as Equation (63) were unsuccessful. Numerical studies using the MSM 
guidance scheme therefore were conducted with the assumption of a contin- 
uously correcting procedure. 

An immediate disadvantage is obvious in the guidance scheme 
represented by Equation (62) . Since values for the converged Riccati 
variables are not stored by the MSM except at the initial point .these 
variables must again be integrated forward from this initial time re- 
gardless of when the perturbation 6x(t) occurs. A more detailed dis- 
cussion of this problem is contained in the next chapter. 



CHAPTER 5 


DISCUSSION OF NUMERICAL RESULTS 
The modified sweep method algorithm was programmed for the 

l 

\ 

LMIVAC 1108 digital computer at the Manned Spacecraft Center in Houston, 
Texas. The integration schemes used follow. 

5 . 1 Numerical Integration Routines 

Fixed Step-Size Integration . Fixed step-size integrations were 
carried out using an Adams-Bashforth predictor-corrector procedure with a 
Runge— Kutta starter (Eastman and Fowler 15 ). The Adams predictor had a 
discretization error of o(h 5 ), and the Bashforth corrector discretization 
error was of o(h 6 ); h is the step-size. The Runge-Kutta starter had a 
discretization error of o(h 5 ). Partial double-precision arithmetic was 
used as follows: the values of the dependent variables were carried in 

full double precision, but the derivatives were evaluated and stored as 
single-precision numbers. This technique minimized the effect of round- 
off error. 

Variable Step-Size Integration . Variable step-size integrations 
were carried out using a predictor-corrector-starter procedure as mentioned 
above. -However, the discretization error in all cases was of o(h 5 ). 

These integrations were carried out in full double precision (Schwauscn 33 ) . 

/ 

5.2 A Brachistochrone Problem 

To compare Modified Sweep Method converged results with known 
analytical solutions for a problem of sufficient complexity, a class of 
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free- final- time Brachistochrone problems was posed as follows. 



Minimize the value of the final time t^ for a particle to fall 
along a frictionless path in a constant gravitational field from point 1 
to point 2 subject to the constraints 


t-* 


X = 

V cos u 

• 

y 

V sin u 

rt 

O 

II 

0.0 

y(t o ) = 

1.0 


where 


and 


V 


a 


/ 2g(y - a) 

v o 2 

y o"2i 


me variational Hamiltonian is H = V(A cos u + A sin u). 

x y 


Two differ- 
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ent cases were scl / - 

r E '1r the terminal constraint vector : 


f /{ h £ ) -5.0 = 0 


and 


M 2 a f'(t-) - 5.0 
I ! £ 

!/ (t f ) - 8.0 


V 


The Modified Swee- 


1 l Ionian and its Partial Derivatives 


Using the . T 

duality condition H - 0 and the strengthened 

Legendre-Clebsch c; T 

Mf, cn that <Su H <5u > 0 , the control vector can 
uu 

be eliminated to 

Hie following Euler-Lagrange equations 


R 

. i 

. i V 

y\ I 

i j 

i 

fx 1 

<* ^ 

. / 

0 

K 1 

.8. 




~ T 
= H 

X 


= — H 


x 


where 


A i _ 

Furthermore, the se - . 


+ X 2 


• 


partial derivatives required by the perturbation 


equations are 


H — - i 

xx I , 

< 1 r * 


0 

Asi 

V 3 



4 ? 


Ax 


0 


0 



_ £. ( JL 

V ^ A 


) 

) 



The and M| cases were computed using the fixed step-si 2 e 

integrator mentioned on page 44 . For comparison purposes , these two 
cases were solved lasing the Method of Perturbation Functions program, 

MPF (see Lewallen^). A step— size h = 0.01 second was used in all 
cases, with the initially assumed values of the unknown Lagrange multi- 
pliers and final time as follows : 


A° - -0.236§ sec/ft 
= -0,6095 sec/ft 
• t^ - 0.5410 sec 

—8 

The convergence criterion e was specified as 0.1 x 10 . The correc- 

tion. procedure used in all cases was 25%, 50%, 75%, and 100% from the 
fourth iteration onward. 

Rate of Convs rgs nee, Figures 1 and 2 show plots of the 


terminal constraint norms versus time for the Brachistochrone problem. 
Both the MPF and MSM results are plotted. Figure 1 shows the error 



NORM OF THE TERMINAL CONSTRAINTS 


l.E+01 

1.E4-00 

1.E-01 

1.E-02 

1.E-03 

l.E-04 

1.E-05 

l.E-06 

l.E-07 



NORM OF THE TERMINAL CONSTRAINTS 
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FIGURE 2 


M 2 CASE FOR BRACHISTOCHRONE PR03LEM 
f 
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norm for the case M^ = x(t^) - 5.0. The two methods converge m seven 

iterations; however, the MSM consistently shows a smaller error than 

the MPF for each iteration. Note that the decrease in the error norm 

for the MSM is significantly less than the MPF for the last iteration. 

Figure 2 shows the error norms for the case M^ = (xCtj.) - 5.0 

y(t^) - 8.0J ; the MSM error norm is not always less than, the MPF 

error norm. However, the error difference is never large. The terminal 

stages of iteration reveal the same high rate of convergence as in the 

\ 

Ml case, 
f 

•S. 

The MSM for this problem at its worst took 20% less compu- 
tational time. However, this figure is not considered significant be- 
cause two unrelated computer programs were used. 

Accuracy of Converged Results. The MSM converged solutions 

gave six decimal place agreement for both the Mj, case and the M^ case 
when they were compared to the known analytical solutions. For the M^ 

e o q 

case, the initially guessed multipliers X and X were in error with 

x y 

the converged values 224% and 282%, respectively, the initial guess at 

the final time was 2% in error. For the M| case, the initial guesses 

on the same X° and X° multipliers were 500% and 260%, respectively, 
x y 

In this case, the initial-guess error for the final time was 14%. These 
results are tabulated in’ Table IV. 

- To summarize, the MSM has exhibited rapid terminal convergence 
and reasonable convergence envelopes fcr the free- final- time Brachisto- 
chrone problem. 
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TABLE IV 

MSM INITIAL-GUESS ERROR RESULTS 



M f 

Converged 

Values 

% Initial 
Error 

. 

Converged 

Values 

*1 

% Initial 
Error 

x° 

X 

-0.0689 

244% 


-0.0357 

/ 

500% 

x° * 
y 

-0.1623 

282% 


-0.1726 

260% 


0.5271 

2% 


0.6277 

14% 


The excellent results warranted further applications of the MSM 
to more complex problems x<rhose analytical solutions were not known and 
which were of current interest to the space program. For these reasons, 

an Apollo three-dimensional reentry problem was chosen. 

» • 

5.3 Apollo Three-Dimensional Reentry Problem 

In the time interval t < t < t^ , find the roll angle program 
$(t) which can be used to control an Apollo spacecraft so as to minimize 
the weighted sum of heating and acceleration effects 



/l 2 + D 2 . 

-r 

m 


„ 1 

X p 2 V 3 

O 


]dt 


I 
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where ^ X q is a constant weighting parameter. Here, the first-term m the 
integrand serves to measure the acceleration effects due to aerodynamic 
forces while the second term measures the convective heating experienced 
by the spacecraft. This minimization is to be accomplished subject to the 
differential equations of motion given as follows 


f • \ 

h 


„ . 'l 

V sm y 

♦ 

e 


V cos y cos A/ (R + h) cos A 

* 

A 


V cos y sin A/ (R + h) 

« 

V 


G sin Y _ D 

« 

Y 


(G cos y/V) + (v cos y/ (R + h)) + (R cos 8/V) 

* 

A 

J 


(-V cos Y cos A taii A/ (R + h)) - ((L sin 8/(V cos y)) 


The following initial conditions represent the reentry conditions 

</ 

for a space vehicle on an Apollo- type lunar return mission. 


rt 

o 

s 


400,000 ft 1 


'75.757576 mi 

e(t ) 


0.0° 


0.0 rad 

0 





A(t o ) 


0.0° 


0.0 rad 

v(t o ) 


35,000 fps 


6. $181818 mi/sec 

Y(t o ) 


-6.5° 


-0.11344640 rad 

A(t o ), 


0.0° 


0.0 rad 


where 

G = -y/(R + h) 2 

D = pSV 2 C D /2m 
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and 


terminal 


L = pSV 2 C L /2m 
* 

-eh 

p - p e 
o 

m = spacecraft mass (assumed constant) 

i 

Optimal reentry trajectories were determined for two cases of 
conditions 


Mi 

r 


h(t f ) - h f 

e(t f ) “ ® f 

6(t f ) - A f 
V(t f ) - V f 
Y(t f ) “ Y f 


and 


A(t f ) - A f 
V(t f ) - V £ 
A(t f ) - A f 


Values for the terminal state represent a typical set of conditions at 
drogue parachute deployment for the Apollo space vehicle 


V 


'75,504 ft 1 



24.1° 

Af 


-0.6° 


— 

856 fps 



-44.3° 

Af. 


-29.4° 

W 7 


Numerical values for the Apollo parameters were assumed as 



54 


S = 129.3 ft 2 = 0.46379993E T 05 mi 2 

m - 204,0 slugs 

P o » 0.27E-02 slug/ft 3 - 0.39743447E + 09 slug/mi 3 

| = 0.42E-04 ft" 1 = 0. 2217600E + 00 mi" 1 

s 

~ 1/2 

X = 0.24509804E-07 sec/(slug-ft ) ' 

° 1/2 

= 0.17809708E-05 sec/ (slug-mi) 

y = 0. 14076519E + 17 ft 3 /sec 2 = 0.95629856E + 05 mi 3 /sec 2 

R = 0. 20908800E + 08 ft - 0.39600000E + 04 mi 


Figure 3 shows the essential geometrical relationships between 
the state variables for the three-dimensional Apollo reentry problem. The 
variables chosen to specify the state of the point mass spacecraft were 
h = altitude, 0 = longitude, A = latitude, V = speed, y - angle of at- 
tack, and A = heading angle. The following assumptions have been made: 
the earth is a nonrotating homogeneous sphere with its center fixed in 
interial space. Furthermore, its gravitational potential is characterized 
by an inverse-square law and it possesses an exponential atmosphere. 


The Modified Sweep Method Reentry Hamiltonian . In Appendix C, 
the Apollo reentry optimization problem is restated. The mechanics of 
restructuring the problem Hamiltonian by use of local control optimality 
and the strengthened Legendre-Clebsch condition are shown. The Hamiltonian 
which is optimal with respect to the choice of roll angle 8 is given as 
follows : 
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H * (cSV 2 /c 2 + C 2 / 2 m) + X p 1/2 V 3 + X V sin y + 


L D 


(R + H) 


‘ f CO f -J (X - X sin A) + (X sin A + X )} H 

. ^ C ° S A 2 6 3 s'J (R + h) 2 

X cos y 

’ f X sxn y -f 1 

14 v J- 


_ 1 


~ pSV fx C^V + — A 2 + X? cos 2 y] 

2m [ 4 D cos y 6 5 


The Euler^Lagrange equations for this problem are then generated 
by taking fir^t partial derivatives of H as follows : 

• _ Bit • . 3H 

x i " 3^ and X = -— where i,j = 1, ...» 6. 

I J - j 

These results, along with the second partial derivatives H , H and iL. 

' XX XA A A 

which serve as coefficients for the matrix Riccati equation are also pre- 
sented in Appe U di x C. 

Results . Initial attempts to solve the Apollo three-dimensional 
reentry problem encountered some difficulties when the modified sweep 
method algorithm was used. Using the system of units ft/lb/sec, certain 
elements of tU e Riccati matrix K grew very large at the initial time. 
Because of thu se large values, the matrix F also became very large at 
the initial time. Consequently, when its inverse was used to compute 
initial-time Corrections, they were so small that the initial values al- 
tered only in the seventh decimal place. As a result, the initial tra- 
jectory was essentially duplicated by subsequent sweeps. 
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An attempt to impose arbitrary bounds on these variables was 
tried for several bounding orders of magnitude ranging from 0.5 to 
1 x 10°. In all cases , every element in the Riccati matrix achieved the 
bounding value by the third sweep. An attempt was then made to generate 
more accuracy by cycling through the evaluation-correction procedure (EC)^ 
of the fixed-step-size integrator several times. Values were tried for N 
ranging from 2 through 9. This effort to prevent the Riccati matrix 

from going onto the limiting boundary was unsuccessful. 

A scheme which used a scaled fractional part of the corrections 

6 \ was then attempted, and this did not eliminate the numerical diffi- 
culties with the Riccati matrix*. The vector was then altered wiuh 

respect to size and to choice of terminal state variables, neither of 
which was successful. The system of units then was altered to slug/mi/sec, 
for which the range of Lagrange multiplier magnitudes became smaller. Al— 
tering the unit system was tried after discussion with Williamson , 
whose studies on the same problem with the MPF revealed a correlation 
between the numerical sensitivities of the Lagrange multiplier equations 
and the unit system chosen. The choice of slug/mi/sec achieved a more 
suitable scaling for the magnitudes of the multipliers; however, this did 
not Succeed in eliminating the difficulties with the Riccati matrix. 

A variable step-size integrator routine was then introduced 
which revealed the numerical sensitivity of the Apollo three-dimensional 
reentry problem to the single-step error on the UNIVAC 1108. This sensi- 
tivity was measured by fixing the final time at t^ = 437.263 seconds; 
the initial values for the state and Lagrange multipliers were defined as 
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Y 

o A as shown in Table V. The state and Lagrange multiplier dif— 

o \ 

^ er Sntial equations were then integrated forward from t = 0 seconds to 

t f' Using the terminal values for the state and Lagrange multipliers, 
the 


out . 


process was reinitialized at and a backward integration carried 

The values obtained at t Q using this procedure were then compared 


agreement with the defined values of X and A . Four cases were 

o o 

in which the single step error e was bounded: 


Case X 1.0 x lo“ 10 < e < 1.0 x lO -07 

Case II 1.0 x 10 _12 < e < 1.0 x 10“° 9 

Case III 1.0 x 10 _13 < s < 1.0 x 10“ 10 

Case IV 1.0 x 10 _14 < e < 1.0 x lO** 11 

To 

numerically integrate forward from the" initial to the final time and 
re3 hitialize and integrate backward to reproduce the initial values to 
e3 ‘ 7 ht decimal places, the single-step error e had to be bounded as 

-14 -11 

1 X 10 < £ < 1 X 10 A 


'^Vn the error became less than 1 x 10 the integration, step size 

—11 

doubled for the next step. If the error exceeded 1 x 10 , the 

^-■p size was halved; otherwise, the step size remained unchanged. These 
'"'"rules are summarized in Table V where the bar under the digits denotes 
^ fc Viations from agreement with initially assumed values. 

Figures 4 through 9 give a particular set of numerical results 
rc, t this problem. This set of results was essentially identical for both 

of terminal conditions M* and Mr. . 

i t 
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TABLE V 

SINGLE-STEP ERROR TOLERANCE 
INITIAL VECTOR 


i 


0.71895168 - 02 


0.75757576 + 02 

-0.72100154 + 01 


0 

0.81929784 + 01 

X - 
o 

0 

0.24206058 + 01 


0.68181818 + 01 

0.16453211 + 02 


-0.11344640 + 00 

0.35760665 + 01 


0 


NUMERICAL RESULTS 



Case I 

Case II 

Case III 

Case IV 

EREMIN 

_ 1.0D-10 

1 . 0D-12 

1.0D-13 

1.0D-14 

ESSKAX 

l.OD-07 

1.0D-09 

l.OD-10 

1.0D-11 

X h 

0.71934619-02 

0 . 71895186-02 

0.718952,65-02 

0.71895118-02 

x e 

-0.72100154+01 

-0.72100154+01 

-0.72100154+01 

-0.72100154+01 

X A 

0.81929769+01 

0.81929783+01 

"0.81929784+01 

0.81929784+01 

X V 

0. 24207,833+01 

0.24206019+01 

0.24206013+01 

0.24206059+01 

X Y 

0 . 16458225+02 

0.16453217+02 

0.16453225+02 

0.16453211+02 

*A 

0.35760_746+01 

0,35760661+01 

0.35760665+01 

0.35760665+01 

h o 

0.7576,6004+02 

0.75757126+02 

0.75757596+02 

0.75757571+02 

0o * 

-0.40744522-06’ 

-0.13346897-01 

-0.15685918-08 

-0.14599980-09 


-0.65481932-06 

-0.13262769-08 

-0.18274255-08 

-0.10809805-09 

v o 

0.68185440+01 

0.68181816+01 

0.68181827+01 

0 . 68181819+01 

Yo 

-0.11346347+00 

-0.11344641+00 

-0.11344641+00 

-0.11344640+00 

A 

0.41275475-05 

0.20899670-17 

0.10682438-07 

0.73654477-09 

o 
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Figure 4 shows the altitude and scaled reentry speed histones 
that the Apollo spacecraft would follow during optimal reentry for both 
the and cases* Figure 5 shows the histories for the longi- 

tude a latitude, angle of attack and heading angle state variables during 
these optimal reentries. These two figures are of interest because they 
define the optimal state histories which the Apollo spacecraft should fly 
to achieve the specified terminal conditions while minimizing aerodynamic 

acceleration and convective heating. Figures 6 and 7 show the A, 

h 

and A^ multiplier histories, respectively. These are the Lagrange 
multipliers associated with the rates of change in altitude and latitude. 
They are included here to define the trends to be anticipated for the 
specified set of initial reentry conditions. 

The two Lagrange multipliers which are required to define the 
optimal reentry roll profile are shown in Figure 8. This particular fig- 
ure shows the A and A. histories where X and A. are associated 
y A y A 

with the rates of change in the reentry angle and heading angle, respec- 
tively. Figure 9 shows the reentry history for the payoff function. 

The aerodynamic acceleration and weighted convective heating experienced 
by the spacecraft have been plotted to reveal their individual character- 
istics. A study of this figure shows that two peaks occur in both space- 
craft acceleration and heating during the optimal reentries. The optimal 
reentry roll procedures seem to call for a trade— off philosophy between 
acceleration and heating experienced by the Apollo spacecraft. The high 
peak in convective heating is initially balanced with a smaller acceler- 
ation peak in the vicinity of 100 seconds. This situation is subsequently 
reversed in the vicinity of 400 seconds where the high acceleration peak 
is balanced with the smaller heating peak. 
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. Figure 10 and 11 show the X^ histories for the and M 2 

cases, respectively. This is the Lagrange multiplier which is associated 
with the rate of change of speed for the Apollo spacecraft. These figures 
were included to reveal, for the specified initial reentry conditions, the 
sensitivities of this variable to a change in terminal conditions. It is 
interesting to note that both histories are similar with the major dif- 
ference arising beyond 400 seconds. 

Figures 12 and 13 show the optimal reentry roll programs for 

the and M f 2 cases, respectively. These are the roll profiles that 

an astronaut would have to use during reentry from a lunar mission to 
minimize aerodynamic acceleration and convective heating while satisfying 
the desired terminal constraints on the vehicle state. 

In each case, optimal reentry calls for the spacecraft to com- 
mence the reentry maneuver with the lift vector pointed almost straight 
downward. The spacecraft is then ..quickly rolled such that the lift vector 
is pointed almost straight up after 90 seconds. A slower downward roll 
of the lift vector is then initiated so that this vector is at a value of 
169 degrees by 350 seconds. The lift vector is subsequently rolled up- 
ward to approximately 15 degrees by 410 seconds at which time terminal 
downward roll procedures differ depending on the specified values for the 
final vehicle state. Specifying two less conditions on the final state of 
the Apollo spacecraft calls for less terminal roll of the lift vector. 

An attempt was made then to obtain a precise evaluation of the 
MSM computational characteristics. For comparison purposes, the 
three-dimensional Apollo reentry M 2 case was chosen. The MSM program 
was then altered to assume computational characteristics similar to 



500 


TIME, SEC 


FIGURE 10 

\ v VERSUS TIME FOR M J CASE 




70 


DEG 



TIME, SEC 


FIGURE 12 

-CONTROL U VERSUS TIME FOR THE M? CASE 

t 


n!H»:3Wy»3»2S 



72 


<% q i 

Williamson's -3 MPF program. Values for the initial Lagrange multipliers 
and final time were specified in varying degrees of accuracy ranging from 
four to eight decimal places. Both programs were run on the University 
of Texas CDC 6600 digital computer using single-precision arithmetic 
everywhere except in the variable step-size numerical integrators where 
partial double-precision was used. Each integrator had a single-step 
error tolerance of l.OE-lO^ e< 1.0E-08. The correction procedure re- 
quired correcting 100% of the terminal error after each iteration. 

Results of this comparison study are summarized m Table VI. The 


TABLE VI 
MPF/MSM 

Convergence Characteristics 
for M| Case of the Apollo Reentry Problem 


Significant Time to Converge % More Time Number of 

Digits for CDC 6600 Required by Corrections 

t^ (Seconds) MSM Required 

MPF MSM MPF MSM 


8 

33 

49 

49% • 

1 

1 

6 

66 

106 

60% 

2 

2 

4 

165 

207 

s 

25% 

5 

4 
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MPF .program required less time to converge in each case. However, with 
decreasing accuracy in the initial guesses for A q and t^, the MSM 
revealed a tendency toward fewer required corrections and more competitive 
time-to-converge. Terminal error norms of both the MPF and MSM pro- 
grams for the case of four significant digits in A^ and t^ are shown 
by Figure 14 . 

No direct computational comparisons with the SSM were avail- 
able. McGregor 22 used the SSM to converge the three-dimensional 
Apollo reentry problem for the case of terminally specified values for 
6, A, and V. This particular case was converged using a fixed step-size 
integration routine. However, as was discussed previously, the MSM re- 
quired a variable step-size numerical integration scheme to preserve the 
numerical integrity of the state and Lagrange multiplier equations. In 
addition, the MSM failed to converge this particular case of the three- 
dimensional Apollo reentry problem. This failure is currently under study 
by this author. Numerical comparison between the SSM and MPF for this 
particular case of the Apollo reentry problem can be found in the study by 

O 

Tapley, et al .^7 

5.4 MSM Guidance Results 

The MSM guidance scheme was implemented for the three-dimen- 
sional Apollo reentry M 2 case. A 5% perturbation in altitude, speed, 
and angle of attack was initiated at t = 0 seconds to study initial 
reentry condition perturbation effects. A similar perturbation was initi- 
ated at t = 75 seconds to correspond to initial peaks in spacecraft 
heating and acceleration. In either case, it was assumed desirable to 
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correct so as to satisfy the originally specified terminal conditions. The 
guidance scheme then assumed that 

y (t) = z(t) = 0 (t Q < t < t f ) (64) 


Control corrections reduced to 

$u = ^U x + U x (K + D» u + A* 21 ))6 x (65) 

Numerical results revealed that the KSM guidance scheme failed 
to satisfy desired terminal conditions for specified vehicle state pertur- 
bations. Investigation revealed that numerical instabilities arising 
from attempts to forward-integrate the matrix Riccati equation were re- 
sponsible for compromising effective terminal guidance. Further study 
to suppress these instabilities is needed to achieve an effective MSM 
guidance scheme. 



CHAPTER 6 


CONCLUSIONS AND RECOMMENDATIONS 

A new second-order variational method (the Modified Sweep 
‘ Method) was developed for solving the two-point boundary value problem 
of trajectory optimization. If differs from the original Successive 
Sweep Method in that the iteration process is now associated with modi- 
fying the initial values of the Lagrange multipliers instead of the 
control function over the time interval of interest. This approach re- 
quires considerably less computer storage and yields the Euler ian control. 
The new method was tested successfully on several classes of problems. 

The following conclusions were reached about the Modified Sweep Method: 

CONCLUSIONS 

1. The method has appeal for problems in which knowledge of 
the Eulerian control is critical. The MSM yields the Eulerian control 
over the entire time interval of interest upon convergence. 

2. Significantly less computer storage than the SSM was re- 
quired. Only 3n + 2(q+l) quantities were required by the MSM algorithm 
to compute the desired corrections from one iteration to the next. This 

is a desirable characteristic for large r-c. miens lonal problems and small- 
storage computers. 

3. The numerical integration of a least n+q less variables 
than the SSM is required. Tnis feature is desirable because less com- 
putation time is required. 
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4. Rapid terminal convergence characteristics of second-order 
numerical optimization methods is retained by the MSM. 

5. The conjugate-point test feature contained in the original 
SSM is also retained by the MSM. 


1 

6. A numerical comparison of the MSM with the MPF for the 
class of free final-time Brachistochrone problems revealed that the MSM 
possesses acceptable convergence envelopes and competitive time-to- 
converge features. 


RECOMMENDATIONS 


1. The basic nature of the generalized Riccati transformation 
technique for solving the linear two-point boundary value problem of con- 
trol optimization should be studied. It is possible that other equivalent 
combinations might possess a better structure for solving the two-point 
boundary value problem than the combination presently being used. 

2. Sensitivity of the MSM algorithm to classes of problems 
should be determined. This recommendation is made because of the method's 
failure to converge the three-dimensional Apollo reentry problem when ter- 
minal state values are specified for longitude, latitude and speed. 

3. The correction procedure used with the MSM should be opti- 
mized such that the largest allowable correction is always attempted during 

a given iteration. This should be accomplished for the obvious reason of 

» 

reducing computational costs by requiring fewer iterations. 

4. Properties of the other Riccati transformation variables 
and their relations to the reference trajectory should be studied. Cur- 
rently, only information about the Jacob i-Mayer conjugate point condition 
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and the abnormality condition is being extracted. This information is 
contained m only two matrices of the many used by the Riccati transfor 
mation technique. 

5. The MSM algorithm should be extended to treat state as 
well as control inequality constraints. The need for this extension is 
obvious since most practical problems are subject to such constraints. 
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APPENDIX A 


THE INHOMOGENEOUS SET OF PERTURBATION EQUATIONS 

Let x(t), u(t), and X(t) be functions associated with an 
extreme trajectory for the functional to be optimized. With the assump- 
tions made in Necessary Conditions , page 12 , the following Euler-Lagrange 
equations are necessarily satisfied: 

x *= H^(x,u,X,t) = f(x,u,t) (A.l) 

X = -H^(x,u,X,t) (A. 2) 

0 = H^(x,u,X,t) - (A. 3) 

where ( ) indicates that the variables are to be evaluated on the ex- 
treme trajectory. 

1 

Now assume a nearby trajectory characterized by the 2n+m func- 
tions x = x + Sx , u = u + Su , and X = X + <$X. Substituting x, u, 
and X into Equations (A.l) through (A. 3) and expanding into a Taylor 
Series to first order about this nearby trajectory, the following equa- 
tions are obtained 

6x — H. & + H, 6u (A. 4) 

Xx Xu 

SX = -H <Sx - H 6u - H .6A (A. 5) 

xx xu xX 

6H T = H 6x + H <5u + H <5X (A. 6) 

u ux uu uX 

where the partial derivatives of the Hamiltonian H are evaluated along 
the nearby trajectory. 
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Making the assumption that the (mxm) matrix H ^ is non- 
singular, the control corrections can be obtained from Equation (A.(>) as 

6u = H -1 [sh T - H 6x - H ,5X) - (A. 7 

uu ^ u ux u X f 

Using this expression to eliminate Su from Equations (A. 4) and (A. 5) 
then gives 


ft 

6x 

__ 

f 

A 

B 


6x 


r \ 

V 




T 





6A 

V J 


-C 

-A 1 




-w 
^ / 


(A. 8) 


where 


A = -H, H _1 E + H. 

AU UU UX AX 


B = 


-H H _1 H , 
Xu uu uX 


C = -R H _1 H + H 

XU uu UX XX 


V - 


Hx u H "6H 


-1 T 
: 6H 
uu u 


w ~ 


-1 T 
H H <$H 

XU uu u 


Equation (A. 8) represents the inhomogeneous set of linear perturbation 
equations used by second-order variational methods. For computational 
purposes, the following is used: 

6H*(t) = -e u H^(t) , 0< e u <l (A. 9) 

BOUNDARY CONDITIONS 


Boundary conditions for these equations are derived in Appendix 
B. They are summarized here as 
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fix(t ) = fix = 0 

o o 


(A. 10) 


and 


fiX(t^) = (P ) fix. + M dv + a,.dt c + t, 
£ xx f f x f f £ 1 


where the vectors a and x are also defined in Appendix B. 


(A. 11) 



APPENDIX B 


FIRST-ORDER PERTURBATION OF TERMINAL CONDITIONS 


To allow for changes in the variable final time from iteration 

to iteration, the following linear approximation is used throughout this 

section. For an appropriate variable; e.g. , v^, assume that 

dv,. = <5v-. + v_ dt , . 
f f f f 

The transversality conditions on the terminal values of the La- 
grange multipliers are expressed by the condition 


-E T = 


A first-order perturbation of this condition gives 


(B.l) 


a h - < P =cdf dx f + ^h dV + <P *t>f d V d> f <B - 2) 

Replacing dx^ and dX^ using the linear approximation stated in the 
first paragraph above, and grouping terms yields 


(dp t 

"f - + c Vf d ' 1 + ( dT - 

• * T 

Replacing X by use of X = — and trahsposing the 6X^ term to the 

left gives 


dt- - 6X- 
f f 


(B.3) 


• 6X f - + v v + 


DP 


Dt 


X + H T Idt - dS, 


(B.4) 


The required terminal conditions on the state variables are 


given by 


M f = M(x f ,t f ) = 0 


(B.5) 
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Similarly, a first-order perturbation in the variables yields 

dM_ = M <5 + M_dt. (B.6) 

£ Xf x f £ r 

The transversality condition on the final value of the Hamilton- 
ian is given by 

0 f - + H(x,u,X,t)^ (B.7) 

After using the linear assumptions stated in the first paragraph of this 
section, a first-order perturbation in the variables gives 




85 


dft { = 


a r - H H” 1 (H + H .P )]. 6x. 
u uu v UX uX XX Jf 1 


!~) - HE L H H T ] d\) 
^ Dt J u uu uX xjf 


-2. (p + H) + x T — - + x T H T - H H _1 H a 
(Dt ' t Dt xu uu uX 


'f dt f 


H - x T dZ f + HB 'V. dS 

u uu u r u uu m f 


(B.10) 


where 


* r , , „ t \ 

“ * (5T + H x]t 


Manipulating the coefficient for the dt f term, it is possible to rewrite 
this coefficient as 


[JL + - sj ~ - H H ^ ot, 

(Dt \Dt H J f Dt u uu uX If 


In matrix form. Equations (B.4), (B.6), and (B.10) can be written as 


r " 

6X f 


( p ) f 
XX f 

T 

M 

x f 

a f 


<5x f 


f 

T 1 

dM, 

= 

M 

X f 

0 

- 


dv 

i 

i 

0 

dt! f 

X J 


T a. 

«f + T 2 

•T 

M £ +T 3 

&f + T 4 


dtf 

1 

1 



(B. 11) 
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t = -dE, 
I f 


-fH H _1 (H + H P )} 

[ U UU UX uA XX 


■fH H" 1 H ,M x 1 

[ U UU UA Xjf 


A 

T 4 


T Df -I 

it — + HHl.a 
f Dt u uu uA 


and 


fH H'VV fx T - H H‘V] dE 

{ U UU ujj y u uu uA J t 



APPENDIX C 


A MODIFIED SWEEP METHOD HAMILTONIAN WITH 


FIRST AND SECOND PARTIAL S OF APOLLO 


3-D REENTRY PROBLEM 


Problem Statement 


Minimize the real functional 


J = 


+ x 

n ° J 


subject to 


- V sin y 

V cos y cos A 
(R + h) cos A 


(R + h) 


cos y sin A 


V . D 

• sin y 

(R + h) 2 m 


•_ . y cos y , V , L cos E 

Y — - — ^ + cos Y 3 

(R + h) 2 v (R + h) mV 


A - - 


V cos y cos A sin A L sin E 


(R + h) cos A 


m V cos y 


and satisfying the end conditions 


C q (a constant vector) 


M(x f> t f ) = 0 


a (qxl) vector 
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where 

♦ 

P 

L 


0 


* 


I » sv \ 

I O^D 


The variational Hamiltonian is 


H = 


2m 


P o e"^ h SV 2 / C 2 + C 2 + X o p^ 2 e“^ h/ ^V 3 + X^ sin y 


, V cos y 
+ <R + h) 


— * (A - X. sin 4) + (X, sin A + X.) 

COS A 2 6 3 5 


(R -f- h)^ 


, . . cos y 

Xj* sin y + x 5 “ip- 


* 

1 "3 h OXT 
-r— p e SV 
2m o 


X, CLV - C T U c cos 8 - X c 
4 D L V 5 6 cos 


sin 3 \ ' 
cos y J 


Partials of the Apollo 3-D Hamiltonian With 


Respect to the Roll Angle 


H e 


8H 

98 


_1_ 

2m 


P 




sin 8 


cos g' 
cos y 


= 0 


-This implies that 


X c sin 8 + X 

D 6 


cos 8 
cos y 


0 
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or 


Thus 


We require 


Thus 


tan 8 = — 


X cos Y 
5 


sin £ = 


±A^ + X| cos 2 y 
5 o 


and 


cos $ = 


-\ c cos B 
o 


±Ag + X| cos 2 y 


Sufficiency Condition 


H 


S3 


1 n -8h r 
0 p B SV C T 

2m o L 


X 5 cos 8 


X 6 sin 8 >| 


cos y 


H es =■ 0 


for a local minimum 


, sxn 8 

X_ cos 8 - X, 

5 6 cos y 


< 0 


or, using H g = 0 


cos y 


±A2 + X? cos ^ 


6 5 


< 0 


Requiring " — < y < ~~ yields 
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0 < cos y < 1 


Therefore, the "+ n sign of the radical must be chosen for optimality xn {5. 
Finally, 


sin 3 


OPT 


/l 2 + X 2 cos 2 y 


cos 3 


“X 5 cos y 


OPT 


v'X 2 


+ X 2 cos 2 y 

o b 


The Optimal Hamiltonian With Respect to 3 


Eliminating sin 3 and cos 3 then yields 


H » -£• P e” eh SVVc 2 + C 2 + X p 1/2 e “^ h/2) V 3 
2m o L D oo 


i •> tt _ • , V cos y 


fcos A 


(X„ - X sin A) 


cos A 2 6 


+ (X 3 sin A + A ) 


r L sin Y + X, 1 

(R + h ) 2 l 5 V J 


* 

1 -3h c „ 

' 2S p o e SV 


|X C V + - /x 2 -r X 2 COS 2 Y 

[ ‘t D cos Y 6 5 



91 


First Partials of H With Respect to x and A 


The Euler-Lagrange equations for the MSM are obtained from che 
first partial derivatives of the Hamiltonian H. These equations are given 
below, where 


{ _ /3H \ /3H \ , . . - , 

X i - - Ita.) ' x i“l3X.) and ^ “ 1. ... 6. 

\ x' ' X' 


H‘ - £3. . 

x 3h 


4 ^ e- Bh SV^/SFT3 - i A l/2 e - B Ch/2) v3 


2m 


L D 2 o K o 


V cos y 

(R + h ) 5 


cos A 


cos A 2 


(A - Ac sin A) + (A- sin A + A.) 


2p 


(R + h)' 


A^ sin y + 


A 5 'cos y ) 


■* *£* 4h " 


A C V + -A 2 + A 2 cos 2 y 

4 D cos y 6 5 ' 


H = § = 0 

- X 2 , 96 .. 


H - 
x 3A 

3 


V COS Y 
(R + h) 


cos A 


cos 2 A 


(X„ sin A ~ A ) 


H 


3H 


3V . m 


= — e eh SV/C 2 + C 2 + 3A p 1/2 e & ^ Z \2 
~ L D oo 


+ X 2 sin y + 

(R + h) 


cos A 


(A - A, sin A) 


cos A 2 6 
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+ (X sin A + X,. ) 


Xi- cos y 


(R + h) 2 V 2 


- — p e ® h sf2X, C n V + - - A| + X| cos 2 y 

2m 'o I 4 D o o 

v cos y 


„ 3H , „ V sm Y 

H = r— = X. V cos Y “ VbTT T 

X- 3Y 1 (R + h) 


cos A ,, .. - 

7 O‘o “ *c sxn A) 

cos A 2 o 


+ (X 3 sin A + X 5 ) 


(R + h) 2 


X^ cos y — y“ sin y 


* f 

1 „ 0 "®«Tr 

p <, e SVC L 


X 2 sin y 
6 


v cos 2 y V^X 2 + cos 2 y 
b 


H 


3H 

3A 


V cos y 

(R + h) 


slIt ^ (X 0 - X c sin A) + X o cos A 
cos A v 2 6 3 




V sin y 




V cos y cos A 
(R 4* h) cos A 



_9H_ | _ V cos y sin A 
3X 3 J " (R + h) 




(R + h)‘ 


sm y “ 


_1_ 

2m 


p o e 


& 

-eh 


sv 2 c 


D 
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(R + h) 


(R + h) 


_ COS 

2 V 


1 -gh, M7r , 
2m p o e SVC L 


X 5 cos y 
A| + X 2 cos 2 y 


V cos y cos A 'sin A 
(R + h) cos A 


i _ _4h 5Yc l \ 

2m 9 o“ cos y rf ~ , 2 — 1 


X z + X| cos z Y 


Matrix 


0 0 


0 0 


10 0 


; fc 

C v SVe~ eh 


^X 2 + X 2 cos 2 y J' 


T A 6 COS Y 

\ ^ / o & 


C L SVe 


( X 6 + A 5 C ° s2 


* 

C L SVe*' 0h 


( A 6 + X 5 cos2 ^) 3/2 


A^ cos y 
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H, Matrix 
Ax 


H, = 0 

X i x i 


B, ■= 0 

X l x 2 


H = 0 

X l x 3 


H, = sin y 

X l x 4 


H. = V cos y 

i X 5 


H, = 0 

X 1 X 6 


H 


X 2 x l 


V 

cos 

(R + h) 2 


cos A 
cos A 


H 


X 2 x 2 


= 0 


H 


X 2 x 3 


(R + h) 


cos y cos A sin A 


cos" 


H 


X 2 x 4 


cos y cos A 
(R + h) cos A 
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V sin y cos A 
(R + h) cos A 


V cos V sin A 
(R + h) cos A , 


(R + h)' 


cos y sin A 


(R + h) 


sin A 


(R + h) 


sin Y sin A 


(R + h) 


cos Y cos A 


Vi 


^-3 sin , + ! ^ SV 


H, -0 

X 4 x 2 


0 
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o * 

— SVC 
m D 


(R + h)‘ 


cos y 


V cos y ^ 2y 


(R + h) 2 (R + h) 3 v 


* *o -6h 

+ ^ 2m C L SVe 


X 5 cos y 


'X 2 + X? cos 2 y 

6 5 


cos y y_ 


(R + h) (R + h) 2 V 2 


2m C L^ Bh 


X 5 cos Y 

^X 2 + X 2 cos 2 Y 
6 5 


V sin y ^ y sin 

(R + h) (R + h) 2 V 


^C T SVe-K x 
2m L 5 2 


0 
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where 


A^r sin y 


A 2 + A 2 cos 2 y) 


6 5 


* 

-Bh 


K. 


Vl 


• r , . . . p C T SVe ... 

V cos y cos A sin A + * __o _L & __ 

cos y *A 2 + X 2 cos 2 y 

6 5 


(R + h) 2 cos A 


H. =0 

X 6 X 2 


H 


Vs 


JV cos y cos A 


H 


Va 


(R + h) 2 

cos z A 


cos y cos A sin A 
(R + h) cos A 


p - C_ Se 
o L 


* 

-eh 


2m COS Y Ai + X2""c^2 


6 5 


H. 


Vs 


P * 

V sin y cos A sin A o „ „„ -3h, 
iTiO 2m C L SVe Vj 


H 


v« 


V cos Y sin A sin A 
(R 4- h) cos A 


where 


^X 2 + 2X 2 cos 2 ^ 

c ° s2 Y k +j 5 c ° s2 v ) 3/2 


- s ^ n X 



H Matrix 
xx 


H 


X 1 X 1 


= a 2 — ^ e“^ h SV 2 /c 2 + C 2 -f ~ X p^e^^^Y 3 
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L D 4 o o 


, 2V cos y 

'T 

(R + h ) 3 
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cos A 2 6 3 5 


6y 
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i - . , cos y 
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*■ a 


2m 


* 
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, X C V + — /x 2 + X 2 cos 2 Y 

( 4 D cos Y 6 5 
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X 1 X 2 
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X 1 X 3 


V cos y 
(R + h) 2 


cos A 


(X sin A - X ) 
cos 2 A 2 D 


H 


X 1 X 4 
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e 
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4 D 
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H = 0 

X 2 X 3 
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* X 2 X 6 


H 


X 3 X 1 
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